#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
//#include <magick/ImageMagick.h>
#include <magick/api.h>

#include <sys/shm.h>
#include <unistd.h>
#include <sys/mman.h>


#define PixelIntensityToQuantum(pixel) ((Quantum) \
  (0.299*(pixel)->red+0.587*(pixel)->green+0.114*(pixel)->blue+0.5))


extern int diffuse(float *data, float *phi, int columns, int rows);
extern int initphi(float *phi, int columns, int rows);

Image * openImage(char * filename, ExceptionInfo * exception)
{
  ImageInfo *image_info;
  Image *image;
  image_info=CloneImageInfo((ImageInfo *) NULL);
  strcpy(image_info->filename, filename);
  image=ReadImage(image_info, exception);
  if (exception->severity != UndefinedException)
    CatchException(exception);
  if (image == (Image *) NULL)
    exit(1);

  fprintf(stderr, "\nImage Info:\n"
	  "name: %s\n"
	  "colorspace : %d\n"
	  "type : %d\n"
	  "columns: %ld rows : %ld\n"
	  "depth: %ld colors: %ld\n"
	  ,
	  image_info->filename,
	  image->colorspace, 
	  GetImageType(image, exception),
	  image->columns,
	  image->rows,
	  image->depth,
	  image->colors
	  );

  return image;
}


int process(Image *image, ExceptionInfo *exception)
{
  int i;
  float * d, *phi,pmin, pmax;
  PixelPacket * ipix;
  Image * oimage;
  ImageInfo *image_info;
  unsigned char *oim;
  oim = malloc(sizeof(unsigned char) * image->columns * image->rows);

  d = (float *)malloc(sizeof(float)*image->columns*image->rows);
  phi = (float *)malloc(sizeof(float)*image->columns*image->rows);
  
  ipix = AcquireImagePixels(image, 0, 0,
			    image->columns ,image->rows,
			    exception);

  for(i=0;i<image->columns*image->rows;i++)
    {
      d[i] = ScaleQuantumToChar(PixelIntensityToQuantum(ipix+i));
    }
  initphi(phi, image->columns, image->rows);
  diffuse(d, phi, image->columns, image->rows);

  pmin=pmax = phi[0];
  for(i=0;i<image->columns*image->rows;i++)
    {
      if(pmin>phi[i]) pmin=phi[i];
      if(pmax<phi[i]) pmax=phi[i];
    }
  for(i=0;i<image->columns*image->rows;i++)
    {
      oim[i] = 255*(phi[i]-pmin)/(pmax-pmin);
      //>0?255:0;
    }

  fprintf(stderr,"Here\n");
  oimage = ConstituteImage(image->columns, image->rows,
			   "I",CharPixel,oim,exception);
  oimage->filename[0]='p';
  oimage->filename[1]='0';
  oimage->filename[2]='0';
  strcpy(oimage->filename+3, image->filename);
  fprintf(stderr, "Writeing file: %s\n", oimage->filename);
  image_info=CloneImageInfo((ImageInfo *) NULL);
  WriteImage(image_info, oimage);
  image_info=DestroyImageInfo(image_info);
  free(oim);
  free(d);
  free(phi);
  return 0;
}


int main(int argc, char **argv)
{
  ExceptionInfo
    exception;
  Image *image;

  InitializeMagick(*argv);
  GetExceptionInfo(&exception);
  image = openImage(argv[1], &exception);
  process(image, &exception);

  DestroyExceptionInfo(&exception);
  DestroyMagick();
  return 0;
}
